rm(list = ls())
require(stats4)

data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)



multinom=function(vec){
return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
}

dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)

#some descriptive stats and graphs
expframe=rbind(dataframe21,dataframe22,dataframe23)
expframe=subset(expframe,uid>800)
expframe$correct <-(expframe$state==4&(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9))|(expframe$state==3&(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5))
expframe$incentive <- (expframe$qid==5)*5+(expframe$qid==1)*40+(expframe$qid==6)*70+(expframe$qid==7)*95
accuracy=c(sum((expframe$incentive==5)*expframe$correct)/sum(expframe$incentive==5),sum((expframe$incentive==40)*expframe$correct)/sum(expframe$incentive==40),sum((expframe$incentive==70)*expframe$correct)/sum(expframe$incentive==70),sum((expframe$incentive==95)*expframe$correct)/sum(expframe$incentive==95))
incentivec=c(5,40,70,95)
incentive=c(5,40,70,95)
uidlist=unique(expframe$uid)

#check N
#**
length(uidlist)
#**

indioutframe=data.frame(uid=vector('numeric'),likelin=vector('numeric'),likepow=vector('numeric'),likegen=vector('numeric'),likeinat=vector('numeric'))
for(indii in 1:length(uidlist)){
indiframe=subset(expframe,uid==uidlist[indii])

totcor=c(sum(indiframe$correct*(indiframe$incentive==5)),sum(indiframe$correct*(indiframe$incentive==40)),sum(indiframe$correct*(indiframe$incentive==70)),sum(indiframe$correct*(indiframe$incentive==95)))
totinc=c(sum((1-indiframe$correct)*(indiframe$incentive==5)),sum((1-indiframe$correct)*(indiframe$incentive==40)),sum((1-indiframe$correct)*(indiframe$incentive==70)),sum((1-indiframe$correct)*(indiframe$incentive==95)))

multinom_const=(multinom(c(totcor[1],totinc[1])))+(multinom(c(totcor[2],totinc[2])))+(multinom(c(totcor[3],totinc[3])))+(multinom(c(totcor[4],totinc[4])))

#Begin Estimation Functions
#########################################
#Shannon

acculin=function(kappa){
accuvec=vector('numeric')
for(i in 1:4){
cost=function(accu){
return(kappa*(accu*log(accu)+(1-accu)*log(1-accu)))
}
negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}
optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

likelihoodlin=function(kappa){
accuvec=acculin(kappa)
like=multinom_const+sum(totcor*log(accuvec)+totinc*log(1-accuvec))
return(-like)
}

if(indii %in% c(11,13,14)){
modellin=mle(likelihoodlin,start=list(kappa=0.55), method="L-BFGS-B",lower=c(0.5), upper=c(1000))
}else{
if(indii %in% c(19,31,36,45,49,52)){
modellin=mle(likelihoodlin,start=list(kappa=1.2), method="L-BFGS-B",lower=c(1.1), upper=c(1000))
}else{
modellin=mle(likelihoodlin,start=list(kappa=20), method="L-BFGS-B",lower=c(0.5), upper=c(1000))
}
}
lllin=logLik(modellin)[1]


################################
#General Entropy Model

accugen=function(kappa,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
if(rho==1){
totcost=kappa*(accu*log(accu)+(1-accu)*log(1-accu))
}
if(rho==2){
totcost=kappa*(-1/2*(log(accu)+log(1-accu))-2*log(2))
}
if(rho!=1&rho!=2){
totcost=kappa*(2^(1-rho)/((rho-1)*(rho-2))*(accu*accu^(1-rho)+(1-accu)*(1-accu)^(1-rho))-1/((rho-1)*(rho-2))-log(2))
}
return(totcost)
}

negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

likelihoodgen=function(kappa,rho){
accuvec=accugen(kappa,rho)
like=multinom_const+sum(totcor*log(accuvec)+totinc*log(1-accuvec))
return(-like)
}

if(indii %in% c(13,14,45)){
modelgen=mle(likelihoodgen,start=list(kappa=0.55,rho=1), method="L-BFGS-B",lower=c(0.001,0.001),upper=c(1000,100))
}else{
if(indii %in% c(19,31,36,49,52)){
modelgen=mle(likelihoodgen,start=list(kappa=1.2,rho=1),fixed=list(rho=1), method="L-BFGS-B",lower=c(1.1,0.001),upper=c(1000,100))
}else{
modelgen=mle(likelihoodgen,start=list(kappa=0.21,rho=14), method="L-BFGS-B",lower=c(0.001,0.001),upper=c(1000,100))
}
}


llgen=logLik(modelgen)[1]


####################################
#Power mutual information

accupow=function(kappa,rho){
accuvec=vector('numeric')
for(i in 1:4){

cost=function(accu){
totcost=(accu*log(accu)+(1-accu)*log(1-accu)-2*0.5*log(0.5))
if(rho>0){
return(totcost^rho)
}else{
return(-1*totcost^rho)
}
}

negpayoff=function(accu){
return(-1*(0.4*incentivec[i]*accu-kappa*cost(accu)))
}

optimout=optimize(negpayoff, maximum=FALSE,upper=0.99, lower=0.51)
accuvec[i]=optimout$minimum
}
return(accuvec)
}

likelihoodpow=function(kappa,rho){
accuvec=accupow(kappa,rho)
like=multinom_const+sum(totcor*log(accuvec)+totinc*log(1-accuvec))
return(-like)
}

if(indii %in% c(1,2,4,6,12,22,24,27,44)){
modelpow=mle(likelihoodpow,start=list(kappa=7000,rho=4),fixed=list(rho=4), method="L-BFGS-B",lower=c(0.001,-100),upper=c(10000,100))
}else{
if(indii %in% c(8,15,32,39,48,50)){
modelpow=mle(likelihoodpow,start=list(kappa=90,rho=1), method="L-BFGS-B",lower=c(0.001,-100),upper=c(10000,100))
}else{
if(indii %in% c(14)){
modelpow=mle(likelihoodpow,start=list(kappa=0.5,rho=1), method="L-BFGS-B",lower=c(0.001,-100),upper=c(10000,100))
}else{
if(indii %in% c(19,29,31,36,41,42,49,52)){
modelpow=mle(likelihoodpow,start=list(kappa=0.45,rho=1),fixed=list(rho=1), method="L-BFGS-B",lower=c(0.45,-100),upper=c(10000,100))
}else{
modelpow=mle(likelihoodpow,start=list(kappa=7000,rho=4), method="L-BFGS-B",lower=c(0.001,-100),upper=c(10000,100))
}
}
}
}

llpow=logLik(modelpow)[1]

####################################
#Innattentive
llinat=multinom_const+sum(totcor*log(0.5)+totinc*log(0.5))

addframe=data.frame(uid=uidlist[indii],likelin=lllin,likepow=llpow,likegen=llgen,likeinat=llinat)

indioutframe=rbind(indioutframe,addframe)


print(indii)
}

indioutframe$AIClin=indioutframe$likelin*-2+2
indioutframe$AICpow=indioutframe$likepow*-2+4
indioutframe$AICgen=indioutframe$likegen*-2+4
indioutframe$AICinat=indioutframe$likeinat*-2

indioutframe$minAIC=0
for(i in 1:length(indioutframe$AIClin)){
indioutframe$minAIC[i]=min(c(indioutframe$AIClin[i],indioutframe$AICpow[i],indioutframe$AICgen[i],indioutframe$AICinat[i]))
}

indioutframe$minAIClin=(indioutframe$AIClin==indioutframe$minAIC)
indioutframe$minAICpow=(indioutframe$AICpow==indioutframe$minAIC)
indioutframe$minAICgen=(indioutframe$AICgen==indioutframe$minAIC)
indioutframe$minAICinat=(indioutframe$AICinat==indioutframe$minAIC)

sum(indioutframe$minAIClin)/52
sum(indioutframe$minAICpow)/52
sum(indioutframe$minAICgen)/52
sum(indioutframe$minAICinat)/52

write.table(indioutframe, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/experiment_2_individuals.csv", sep = ",", col.names = T, append = F)


